Chain motion and viscoelasticity in highly entangled solutions of semiflexible rods 
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Brownian dynamics simulations are used to study highly entangled solutions of semiflexible polymers. Bend- 
ing fluctuations of semiflexible rods are signficantly affected by entanglement only above a concentration c**, 
where c** ~ 10 3 L~ 3 for chains of similar length L and persistence length. For c > c**, the tube radius R e 
approaches a dependence R e oc c _3//5 , and the linear viscoelastic response develops an elastic contribution that 
is absent for c < c**. Experiments on isotropic solutions of .F-actin span concentrations near c" for which the 
predicted asymptotic scaling of the plateau modulus G oc c 7// ' 5 is not yet valid. 



Solutions of long polymers become entangled when the con- 
centration or chain length exceeds a threshhold. The nature 
of "entanglement" is obviously different, however, for random 
walks, rigid rods, and semiflexible threads. It has been pro- 
posed that solutions of semiflexible rods, of length L less than 
or equal to their persistence length L p , may exhibit two differ- 
ent levels of entanglement, in different concentration regimes 
IH S S 0] - a loosely-entangled regime, in which only ro- 
tations and transverse translations are hindered by collisions, 
and a tightly-entangled regime, in which transverse shape fluc- 
tuations are also strongly affected. The crossover between 
these two regimes is expected to be associated with a qual- 
itative change in viscoelastic properties, due to the inability 
of a tightly-entangled solution to rapidly relax stress arising 
from transverse chain deformations. Clear experimental evi- 
dence of "tight" entanglement has been obtained only for so- 
lutions of very long actin protein filaments (F-actin), of length 
L <~ L p 10/im and diameter d ~ 10 nm. The evidence 
comes both from visualization of flourescently labelled chains 
JHQ] and from rheological measurements IE & OH ■ It re- 
mains unclear, however, whether bending fluctuations are ever 
significantly hindered in isotropic solutions of any of a variety 
of other well-studied model systems of semiflexible rods with 



a modified tube model (2J, [3l [13|] in which each chain under- 
goes reptation in a narrow wormlike tube. A scaling argu- 
ment due to Odijk and Semenov (111 [5L 13] predicts a tube radius 

R e oc L p (pLl)~ 3 / 5 for c > c** 

Our simulations use a novel algorithm that was designed to 
allow simulation of Brownian dynamics of arbitrarily thin but 
uncrossable wormlike threads. Each polymer is represented as 
a discretized chain of N inextensible rods and N + 1 beads. A 
periodic cubic simulation cell is initially populated with a ther- 
mally equilibrated solution of wormlike chains, by a Monte 
Carlo growth algorithm. At each step of our dynamical simu- 
lation, a trial move is generated for a randomly chosen chain 
by taking one time step of the Brownian dynamics (BP) algo- 



rithm used in previous work on dilute solutions 1141 115L [1 



A trial move is rejected, however, if it would cause the cho- 
sen chain to cut through any other. Whether or not a move is 
accepted, another chain is then chosen at random, and the pro- 
cess is repeated. In this work, we use an algorithm for chains 
with anisotropic friction, with Cll/Ci = 1/2 EH* where £|| 
and £j_ are longitudinal and transverse friction coefficients, re- 
spectively. Details of the algorithm are presented elsewhere 



L ~ L p 11 111 , for which the average chain lengths and aspect 



ratios are all much smaller than those obtainable with F-actin. 
Simulations offer a potentially important complement to the 
experimental study of these systems, which provide access to 
different information and are subject to different difficulties 
than those encountered in experiments. 

Consider a solution of thin semiflexible rods, each of con- 
tour length L and persistence length L p , with L < L p . Let 
c be the number density of polymers and p = cL be the con- 
tour length per volume. Simple geometrical arguments sug- 
gest the following sequence of concentration regimes J3l : At 
dilute concentrations c < c*, where c* oc L~ 3 , chain mo- 
tion is essentially unhindered. In the loosely-entangled regime, 
c* -C c <C c**, rotations and transverse rigid body transla- 
tions are strongly hindered, but tranvserse bending fluctuations 
are not. In this regime, each chain is tr app ed in a cylindrical 
cage or tube of radius R e ~ l/(cL 2 ) 11211 . Above a thresh- 
hold c** ~ ^jLp/L c*, this cage become narrow enough 
to also hinder thermal bending fluctuations (ll 0, S 01- At 
concentrations c ^> c**, chain motion can be described by 



|17t liSD- Here, we present results for entangled solutions of 
chains with L/L p = 0.25 - 2.0 and N = 10 - 40 rods at 
concentrations cL 3 = — 4000. A 1 mg/ml solution of (hy- 
pothetically) monodisperse F-actin filaments with L = 8fim 
would have cL 3 ~ 2500. 

To characterize the effect of entanglement upon bending 
fluctuations, we have calculated two measures of the trans- 
verse mean-squared displacement (MSD) vs. time for the mid- 
dle bead of a polymer. The quantity (Ad 2 (t)), shown in the 
main plot in Figure [T] is the variance of the distance Ad(t) 
between the middle bead at time t and the closest point on the 
contour of the same chain at an earlier time t = 0. The in- 
set shows (Ar^ ± (t)) = (|r m)J _(i) - r m , ± (0)| 2 ), in which 
Ym,±(t) is the transverse component (transverse to the local 
chain tangent) of the displacement of the middle bead from the 
chain's center of mass. The quantity (Arf n j_(f)) is not sen- 
sitive to center-of-mass diffusion, but only to displacements 
arising from bending fluctuations, and so approaches a finite 
value at long times. At early times, both measures of trans- 
verse MSD increase as t 3 / 4 , as predicted Oj]. With increasing 
concentration, both quantities become suppressed over a range 
of intermediate times, indicating the formation of a tube. 
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FIG. 1: {Ad 2 (t)}/L 2 vs t/r° od for chains with L = L p , where 
r® od = C,±L 3 I (72kBT) is the rod rotation time in dilute soluton. 
Numbers near curves are values of cL 3 (or 'Dilute' for c = 0), while 
numbers in parentheses indicate N. The black dashed line is the pre- 
dicted asymptote at early times, for which (Ad 2 (t)) oc t 3//4 . Inset: 
(Ar 2 mt± (t)} vs t/r° od for L = L p and cL 3 = 'Dilute', 100, 250, 
500, 1000, and 2000. The red dashed line is the result of a slithering- 
snake simulation of pure reptation. I"l8ll 

If each chain were confined to a tube of well-defined radius 
R e over a wide range of intermediate times, (Ad 2 (t)} would 
develop a plateau, with a plat eau value (Ad 2 (t)) ~ iR% Here, 
R 2 is defined, as in Ref. 12011 . as the variance of the transverse 
displacement of the chain from the "center" of the tube (i.e., 
the average chain contour) in either of two transverse direc- 
tions. A plateau could appear in (Ad 2 (t)) even in the rigid 
rod limit, however, due to suppression of transvserse center-of- 
mass motion. The suppression of (Ar^ ± (t)) at intermediate 
times, however, is evidence of hindered bending motion, and 
thus of tight entanglement. In fact, we never observe a clean 
plateau in either quantity. Instead, we see a crossover from 
f3/4 g rowm a ); sm all t to a much slower growth at interme- 
diate times, which becomes flatter with increasing concentra- 
tion and/or chain length (i.e., increasing cL 3 ), with a crossover 
time r e that decreases with increasing c. The suppression in 
(Ar^ j_(t)) is signficant only for cL 3 > 500, suggesting a 
crossover c** ~ 500/L 3 for L = L p . 

For cL 3 = 1000, our results for (Arf n ± (t)) include both a 
plateau at intermediate times and an upturn at the end of this 
plateau. This upturn is mimicked very accurately by the re- 
sults of a separate slithering-snake simulation of pure repta- 
tion of a wormlike chain (the red dashed line in the inset) IU8I1 . 
Pure reptation yields a nonzero transverse MSD (Ar 2 n j_(t)) 
at times less than the reptation time r rep = (,\\L 3 / (ir^ksT) 
because reptation occurs along a curved tube. (Ad 2 (t)) is de- 
fined so as not to be affected by pure reptation, and shows a 
slightly broader plateau than (Ar^ n ± (t)). 

To quantify r e and R e , we have collapsed our data for 
(Ad 2 (t)) in a manner that assumes the existence of a scaling 




FIG. 2: Collapse of (Ad 2 (t)} data. Curves labelled "I" represent 
collapse of data for N = 20, cL 3 = 250, 500, 1000, and L/L p = 
0.25, 0.5, 1.0, 2.0. Curves labelled "II" represent collapse of data for 
N = 40, cL 3 = 1000, 2000, 4000, and L/L p = 0.25, 0.5, 1.0, 2.0. 



relationship (Ad 2 (t)) = AR 2 e f(t/T e ). That is, we have cho- 
sen values for R e and an entanglement time r e for each set of 
parameters so as to collapse the data for many different val- 
ues of L/Lp and cL 3 onto a master curve of (Ad 2 (t)) / (2R e ) 2 
vs. t/r e . The resulting collapse is shown in Figure [2] We 
display separate master curves for chains with N = 20 and 
N = 40 because early time behavior is noticeably different for 
discrete chains with different numbers of rods. The collapse 
is excellent for solutions with cL 3 > 1000. The horizontal 
dashed lines with (Ad 2 (t))/AR 2 = 1 represent an assumed 
long time asymptote for hypothetical systems of much longer 
chains, from which we have extracted estimates of R e . 

Figure [3] shows resulting values of the dimensionless tube 
radius R e /L p vs. dimensionless concentration pL 2 for systems 
with L/Lp = 0.25 — 2.0. Dimensional analysis requires that 
the ratio R e /L p be a function R e /L p = f(pL 2 ,L/L p ) of 
dimensionless length L/L p and dimensionless concentration 
pL 2 alone. In the tightly-entangled regime, however, we ex- 
pect R e to become independent of L, implying that R e /L p 
must approach a function of pL 2 alone for c 3> c**. At high 
concentrations, our results for different values of L/L p do 
indeed approach a common asymptote, which is furthermore 
very accurately described by the predicted relation R e /L p = 
a(pLp) -3 / 5 , with a = 0.95 (dashed black line). For each 
value of L/L p , R e /L p also exhibits small but systematic de- 
viations from this asymptote at lower concentrations. This de- 
viation is seen most clearly in the inset, in which we plot the 
ratio i? e /[0.95Lp(pL2)- 3 / 5 ] vs. c/c**, where we have taken 

c** = 500L P ^ 2 L~'/ 2 . The near collapse of the deviations 
from the asymptote for the stiffest 2 chains {L/L p =0.25 and 
0.5) is consistent with the prediction that c** oc Lp^ 2 L~ 7 / 2 for 
L <C Lp. J3l. The experimental values for R e in F-actin so- 
lutions(crosses) are the flouresence microscopy results of Kas 
et al. JHSh as defined and presented previously in Ref. ll20Tl . 
The crossover from loose- to tight-entanglement is expected 
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FIG. 3: Non-dimensionalized tube radius vs. concentration. Num- 
bers in the legend are values of LJL V . Crosses are fiourescence mi- 
croscopy results for F-actin @, 01, non-dimensionalized by L p = 
17/im. Inset: Deviation R e /R^ from the asymptote vs. c/c**, 
where Ri°^ = 0.95L p (pLp)~ 3//5 corresponds to the dashed line in 
the main plot, and where c** = 500Lp //2 L~ 7/ ' 2 . 
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FIG. 4: Non-dimensionalized sum [G or nt(i) + G CU rv(t)]/(ckBT) of 
the orientational and curvature moduli vs t/r^ od for L/L p = 0.5. 
Numbers in the legend are values of cL 3 . Dashed curves are fits, as 
discussed in the text. Black solid curves in both plots are theoretical 
predictions for dilute solutions fiUl . Inset: Corresponding tension 
stress G tens (t)/(ck B T) for systems with cL 3 = 0, 250, 500, 1000. 



to cause a dramatic change in viscoelastic behavior. Detailed 
theories of linear viscoelasticity have been developed for the 
extreme limits of dilute solutions (c <C c*) Il4i 1 1 511 and of 
very tightly entangled solutions (c 3> c**) 111 311 . Both the- 
ories make use of a formal decomposition of the stress into 
curvature, orientational, and tension contributions 01, and a 
corresponding decomposition of the dynamic modulus G{i) 
(i.e., the response to an infinitesimal step strain) as a sum 
G{t) = G curj (t) + G ornt (t) + Gtens(i)- In both dilute and 
loosely-entangled solution, G curv (i) and G te ns(i) are predicted 
to exhibit power law decays at very early times, but to de- 
cay exponentially at times greater than the relaxation time 
tj_ = P^±L /(ksTLp) of the longest wavelength bending 
mode, where = (4.74) -4 . For c < c**, G(t) is thus domi- 
nated at t > t± by a more slowly decaying orientational modu- 
lus G ornt (t) ~ (3/5)c£: B Te-*/ 7W , where r rod is a rotational 
diffusion time. In loosely-entangled solutions, the only pre- 
dicted effect of entanglement is to increase r ro d, without sig- 
nificantly changing G curv (t) or G t ens{t). The plateau of mag- 
nitude (3/5)cfcsT in G orn t(t), which is present even in dilute 
solution, reflects the free energy cost of partially aligning an 
initially random distribution of rod orientations. The crossover 
to tight entanglement, however, is expected to cause a plateau 
to appear in G curv (i), with a plateau value G curv .o that varies 
as liMEl G curv , oc k B Tp 7 / 5 Lp 1/5 for c » c**. 

We have "measured" G(t) and its components by simulat- 
ing stress relaxation after a rapid, small amplitude uniaxial 
step extension of an initially cubic periodic unit cell. Stress 
is evaluated using the virial tensor, as in previous simulations 
of dilute solutions lfl4l [l5ll . Measurements of G(t) in dilute 
solution by this method agree to within statistical errors with 



those obtained previously II 141. Il5ll from stress fluctuations in 
equilibrium. 

In Figure 21 the main plot shows a non-dimensionalized 
sum [Gomt(i) + G curv (t)]/ (cksT) of the two components 
of G(t) that are predicted and observed to exhibit an elas- 
tic plateau. The inset shows Gtens(i)/ (cfcsT), which, as ex- 
pected lfl3ll . does not exhibit a plateau, and which is found 
to be almost independent of c over this range of parameters. 
The plateau in [G ornt (t) + G curv (i)]/ (cfcsT) becomes signif- 
icantly greater than the limiting value of 3/5 obtained in di- 
lute solution, which arises from G ornt (i) alone, only above 
an apparent crossover concentration of c**L 3 ~ 250 — 500, 
above which G curv (i) also begins to contribute to the observed 
plateau. The terminal relaxation is accessed in our simulations 
only for cL 3 < 500, but the plateau value is always accessible. 

To quantify the plateau modulus, we have fit the sum 
Gcurv(^) + G ornt (t) to a function 

f cfcuTe-*/ 7 -"* + G curv , dil (<) + G cunfi e' t/T0 , (1) 
5 

Here, the first term on the r.h.s. is an expression for G ornt (t), 
where T ro( j is a concentration-dependent r otational diffus ion 
time, and G curv ,dii(£) is the prediction of Shankar et alj for 
Gcurv(i) in dilute solution. The quantity G curv ,o is the con- 
tribution of G curv (t) to the overall plateau modulus, which 
is an adjustable parameter. We have used a time constant 
To = T rep /2 for the relaxation of the curvature plateau. This 
was chosen to fit the observed decay of G cmv (t) alone (not 
shown separately here) at cL 3 — 250 and 500, and is con- 
sistent with a double-reptation model of the relaxation of the 
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FIG. 5: (GoLp)/(pkBT) vs pLp from simulations and experiments 
on entangled F-actin solutions. Simulation results for systems with 
L/Lp — 0.5 are shown for both Go (filled squares) and G CU rv,o (filled 
circles). Experimental results for F-actin solutions of Gardel et al. 
IToll (open circles) and of Hinner et al. @l for filaments of unregu- 
lated length (open triangles) and of average length 16/i?n regulated 
by gelsolin (open squares). Dashed line with a slope of 0.4 is the 
prediction of the binary collision approximation of Ref . [20] . 

curvature plateau. Values of the rotational diffusion time 
Trod(c) were measured in separate equilibrium simulations 
111, which yield r rorf /r r O(J = (2.10,3.06,4.46,6.17) for 
cL 3 = (250, 500, 1000, 2000). The values of G curv , obtained 
by fitting this data depend very little upon our choices for the 
time constants tq and T ro d(c). 

The total plateau modulus Go in G(t) is a sum Go = 
(3/5)cfcsT + G curVj o °f orientational and curvature contribu- 
tions. Figure|5]compares simulation results for Go and G clm n 
to reported values of Go in entangled F-actin solutions rfsL TlQfl - 
The results of Hinner et al. JH were obtained by macroscopic 
rheological measurements, while those of Gardel et al. 11011 
were obtained from two-particle micro-rheology. Our results 
for Go agree well with the values of Hinner et al, and are well 
within the scatter of results reported in the recent literature. A 
fit of our results for Go to a power of c yields Go /c oc c 7 . 

It is clear from the simulation data, however, that the range 
of concentrations accessed in our simulations, and most of that 
studied experimentally, lies within about one decade of the be- 
ginning of a broad crossover to tightly entangled behavior, be- 
low which Gcurv(i) does not contribute to Go- As a result of 
this proximity to c**, G ornt {t) dominates Go over much of this 
range, while the contribution G curv o that is actually predicted 
to vary as G curv .o/c oc c 0A in the limit c ^ c** increases 
much more rapidly from nearly zero. The results suggest that 
the very rough agreement between the predicted asymptotic 
behavior of G curVj o for c ^ c** and measurements of Go in 
F-actin may be largely fortuitious. 

The isotropic-nematic (IN) transition for rodlike polymers 
occurs at a concentration cinL 3 ~ 4L/d. Values of L/d for 



available model systems with L < L p other than F-actin, such 
as Fd virus 12211 (F ~ 0.9^m, Fj, ~ 2fim, d ~ 7 nm) and rod- 
like poly(benzyl glutatmate) 11 If (L p ~ 0.l5(j,m and d ~ 2 
nm) are all at least 10 times smaller than for F-actin, for which 
L/d ~ 10 3 . The IN transition in systems with L/d < 100 oc- 
curs at concentrations cinL 3 < 400, at which c < c**. Our 
rough estimate of c** ~ 500Lp^ 2 L~ 7 / 2 for L < L p implies 
that a clear tightly-entangled isotropic regime for semiflexible 
rods can exist only in systems with L/d > 10 3 . This is con- 
sistent with the fact that a clear rheological signature of tight 
entanglement has been observed only in F-actin solutions. 

Taken as a whole, our results both provide evidence for 
the correctness of a simple scaling theory for the asymp- 
totic dependence of tube radius upon concentration in tightly- 
entangled solutions, and (equally importantly) clarify the lim- 
its of validity that theory, particularly as applied to rheology. 
It appears that bending fluctuations of rods with L ~ L p are 
signficantly hindered by entanglement only under surprisingly 
stringent conditions. 
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